% References in [takato % yamaguchi, 1995]
% this function is used for generating the spatial auto-correlation
% for Zernike coefficients based on von-Karman turbulence
% Written by Lu benchu, Taiyuan University of Technology
% email: benchuul@163.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function takato_E_ajs = ...
    takato_spatial_correlated(j1,j2,s,theta0,D,k0,lambda,Cn2,L)
% j1,j2: single-index for zenike polynomials, based on Noll's indexing strategy
% s,theta0: coordinates with dimenson [D,rad(0*pi~2*pi)]
% D: diameters of telescope [m]
% k0: 1/L0 [1/m] ; L0: Outer scale of atmospheric turbulence
% lambda: wavelength [m]
% Cn2(z): refractive structure constant, 'z' with dimension [km]
% L: distance of propagation [m]
R = D/2; % radius of telescope [m]

[n1,m1] = Noll_j_to_nm(j1);
[n2,m2] = Noll_j_to_nm(j2);
A_para = takato_A();
fjj = takato_fjj();
takato_E_ajs = 2^(14/3)*pi^(8/3)*A_para*R^(5/3)...
    *sqrt((n1+1)*(n2+1))*fjj;

% subfunction for calculating the parameter A in takato's equation
% --------------------------------------------------------------
    function A_para = takato_A()
        % lenbda: wavelength
        % Cn2: refractive constant
        % L: length of propagation
        
        if (strcmp(class(Cn2),'function_handle'))
            %if Cn2(1e16) ~= 0
            %    error('error, Cn2(Inf) unequal to 0');
            %else
            funCn2 = @(x) Cn2(x.^10e-3);
            A_para = 0.00969*(2*pi/lambda)^2*integral(funCn2,0,L);
            %end
        else % in cases where Cn2 is a constant
            A_para = 0.00969*(2*pi/lambda)^2*Cn2*L;
        end
        
        % end of the function
    end

% subfunction for calculating f_jj' in Takato's equation
% -----------------------------------------------------
    function fjj = takato_fjj()
        
        m_swit = min(m1,1) + min(m2,1);
        switch m_swit
            case 0
                I_para = takato_I(0,n1+1,n2+1,2*s,2*pi*R*k0);
                fjj = (-1)^((n1+n2)/2)*I_para;
            case 1
                I_para = takato_I(m1,n1+1,n2+1,2*s,2*pi*R*k0);
                if m2 == 0
                    if mod(j1,2) == 0
                        fjj = (-1)^((n1+n2-m1)/2)*sqrt(2)*cos(m1*theta0).*I_para;
                    else
                        fjj = (-1)^((n1+n2-m1)/2)*sqrt(2)*sin(m1*theta0).*I_para;
                    end
                else % m1 == 0
                    if mod(j1,2) == 0
                        fjj = (-1)^((n1+n2-m2)/2)*sqrt(2)*cos(m2*theta0).*I_para;
                    else
                        fjj = (-1)^((n1+n2-m2)/2)*sqrt(2)*sin(m2*theta0).*I_para;
                    end
                end
            case 2
                I_para1 = takato_I(m1+m2,n1+1,n2+1,2*s,2*pi*R*k0);
                I_para2 = takato_I(abs(m1-m2),n1+1,n2+1,2*s,2*pi*R*k0);
                if mod(j1,2) ~= mod(j2,2) % parity: (j1 ~= j2)
                    [n1,m1,n2,m2] = j_parity(j1,j2,n1,m1,n2,m2);
                    fjj = (-1)^((n1+n2-m1+m2)/2)*sin((m1+m2)*theta0).*I_para1...
                        -(-1)^((n1+n2+2*m1+abs(m1-m2))/2)*sin((m1-m2)*theta0).*I_para2;
                else
                    j_mod = (mod(j1,2)+mod(j2,2))/2;
                    fjj = (-1)^(j_mod)*(-1)^((n1+n2-m1+m2)/2)*cos((m1+m2)*theta0).*I_para1...
                        +(-1)^((n1+n2+2*m1+abs(m1-m2))/2)*cos((m1-m2)*theta0).*I_para2;
                end
                
                % end of the switch
        end
        
        function [n1,m1,n2,m2] = j_parity(j1,j2,n1p,m1p,n2p,m2p)
            if mod(j1,2) > mod(j2,2) % j1 is odd
                n1 = n2p;
                m1 = m2p;
                n2 = n1p;
                m2 = m1p;
            else % j1 is even
                n1 = n1p;
                m1 = m1p;
                n2 = n2p;
                m2 = m2p;
            end
            % end of the function for j_parity
        end
        
        % subsub-function for calculating para I in takato's equation
        % -----------------------------------------------------------
        function I_para = takato_I(k,u,v,alpha,x0)
            
            fun = @(x) (x^(-1).*besselj(k,alpha.*x)...
                .*besselj(u,x).*besselj(v,x))./((x.^2+x0.^2).^(11/6));
            
            %if max(fun(Inf)) ~= 0
            %    error('error, the values of I_para is not convergence');
            %else
            I_para = integral(fun, 0, Inf, 'ArrayValued', true);
            %end
            
            % end of the function for takato_I
        end
        
        % end of the function for takato_fjj
    end

% end of the function
end